Partially integrable dynamics of hierarchical populations of coupled oscillators 
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We consider oscillator ensembles consisting of subpopulations of identical units, with a general 
heterogeneous coupling between subpopulations. Using the Watanabe-Strogatz ansatz we reduce 
the dynamics of the ensemble to a relatively small number of dynamical variables plus constants 
of motion. This reduction is independent of the sizes of subpopulations and remains valid in the 
thermodynamic limits. The theory is applied to the standard Kuramoto model and to the description 
of two interacting subpopulations, where we report a novel, quasiperiodic chimera state. 
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Large populations of coupled oscillators occur in a va- 
riety of applications and models of natural phenomena, 
ranging from collective dynamics of multimode lasers and 
Josephson junction arrays to the pedestrian synchrony 
[H; the analysis of the dynamics of these systems is a 
topic of high interest. Even in the context of the simplest, 
paradigmatic case of globally (all-to-all) connected, sine- 
coupled phase oscillators (the famous Kuramoto model 
and its generalizations), many problems remain yet un- 
solved, especially those related to a heterogeneous cou- 
pling and nontrivial collective dynamics. In this Letter 
we treat an important case of a hierarchically organized 
population. It can be viewed as a (finite or infinite) col- 
lection of interacting subpopulations, each consisting of a 
(finite or infinite) number of identical units; sizes of the 
subpopulations and couplings between them are gener- 
ally different (cf. [1,13] and references therein). Using the 
seminal approach of Watanabe and Strogatz (WS) [3], we 
demonstrate that each subpopulation can be described 
by only three dynamical variables plus constants of mo- 
tion, determined by initial conditions. This partial in- 
tegrability allows us to separate the full dynamics into 
a relatively small number of generally dissipative modes 
(their number is proportional to the number of subpop- 
ulations) with possibly nontrivial behavior, and the in- 
tegrals of motion. In particular, our theory allows us to 
extend some recent results [1, Q and to determine the 
conditions of their validity. 

Our basic model is a generalization of the Kuramoto 
model Q, cf. 0,1]: 



re- written as 



dt 
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where Za is the effective force acting on the oscillators of 
subpopulation a. Here we have introduced the relative 
population sizes = Na/N and the complex mean fields 
for each subpopulation 



Xa + iYa 



n: 



fc=i 
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Note that all oscillators in a subpopulation obey the same 
equation, though generally they have different initial con- 
ditions 0^(0). Thus, we can apply to each subpopulation 
the WS ansatz [4j] that reduces the dynamics of the sub- 
population to that of three variables Pa{t), '^a{t), ^a{t), 
via the transformation [8| 
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containing Na constants tp'^, which are directly deter- 
mined from the initial state (/'^(O) and additionally sat- 
isfy 



M Nb 



dt 
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Here we denote the subpopulations by indices a,b = 
1, . . . , M. Variable (/)^(t) is the phase of oscillator k in 
subpopulation a; k — 1, . . . Na, where Na is the size of 
the subpopulation, and Ua is the natural frequency of its 
oscillators (we remind that all oscillators in a subpopu- 
lation are identical). The total number of oscillators is 
N — J2^a and two constants s, a describe the coupling 
with an arbitrary phase shift, cf. The system can be 



Due to an arbitrary shift of constants with respect to 
^P, only Na — 3 of constants ipk ^^'^ independent. The 
WS method is valid generally, provided the number of 
oscillators in a subpopulation is larger than three, and 
the initial state does not have too many clusters, see [3| 
for a detailed discussion of these conditions and of how 
/Oa(0), '^a{0), $a(0), and V'fc can be computed from 4>liO). 
With account of Eq. ([3]), we write the WS equations for 
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our setup as 
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In order to illustrate the physical meaning of the new 
variables, let us consider how they characterize the distri- 
bution of the phases of a subpopulation. Generally, oscil- 
lators form a bunch, and the amplitude p characterizes its 
width: p = 0, if the distribution is uniform (asynchrony) 
and p = 1, if the distribution shrinks to (5-function (full 
synchrony). Amplitude p is roughly proportional to the 
amplitude of the mean field r (see Eq. (jl])) in the sense 
that p = r = for the full asynchrony and p — r = 1 for 
the full synchrony. For intermediate cases these quanti- 
ties generally differ and coincide only in a special case, 
outlined below. The phase variable <i> characterizes the 
position of the bunch, and is therefore related to the 
phase of the mean field, $ f« O. Another phase variable 
^' describes the motion of individual oscillators with re- 
spect to the bunch (generally, the oscillators can move 
with a velocity different from that of the bunch, see 
for an example of such a dynamics). 

The set of Eqs. is a straightforward generalization 
of the WS equations [1] to the case of M interacting 
subpopulations. For a further analysis, and in particular 
for the consideration of the thermodynamic limit, it is 
convenient to introduce new variables, a phase shift — 
— and a complex bunch amplitude Za = PaG^^" ■ 
Then we can rewrite Eqs. ([7][S]) as 



dza _ . 1 „ zl 



dt 



uja + Im(z*Za) 



(10) 

(11) 



Next, we have to represent the complex force Za (see 
Eq. ([3])) in terms of new variables. For this goal it is 
convenient to rewrite Eq. ([5]) in an equivalent form e"^'' = 
e**(pe** -I- e''''^)/(pe*'/''= + e'*). Substituting this into 
Eq. (HI, we obtain: 



Pae^^^^aiZaXa) = Zala{ZaXa) , 
k—1 " 



(12) 



From Eq. pO|l it follows that the dynamics of the com- 
plex bunch amplitude of a subpopulation Za is deter- 
mined by the force Za, resulting from interaction within 
the subpopulation as well as from interaction with other 
subpopulations. Contributions to Za are proportional to 
the relative weights n^, to the coupling constant £e~*" 
and to the complex mean field 7b (z;,, Cb)zb, which gener- 
ally depends not only on the global variables (^i, and Zf,, 



but also on the constants of motion V'fc- Equations (fTOl ^ 
[TT|) . as well as equivalent equations ^]^, together with 
the definitions (|3ll2p are exact and complete; they show 
that the dynamics of a hierarchical ensemble of oscilla- 
tors can be reduced to 3M ODEs plus N — 3M constants 
of motion. Before proceeding with the analysis of these 
equations and examples, let us discuss how a thermody- 
namic limit N ^ oo can be introduced in this picture. 
There are two main ways of performing this. 

(i) Suppose that the number of subpopulations M re- 
mains finite, but their sizes grow N, Na — > oo in a way 
that Ua — const. In this case only Eq. (fT^ is affected 
and should be now written as an integral 



^aiZaXa) = 
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-aaWdi^. (13) 



Here ctq {^p) is the distribution of the constants of motion 
ip in the subpopulation a, additionally it satisfies (cf. (O) 



<JaWe'^ dip = 



(14) 



In this limit the ensemble is described by a set of 3M 
ODEs, where the right hand sides depend on the variables 
via integrals p^ . The integrals of motion are now the 
functions aa{4')- 

(ii) In another limiting case, we keep the size of each 
subpopulation Na finite but let the number of subpopu- 
lations grow M —>■ oo. Considering indices a,b as contin- 
uous variables, we write instead of Eq. ([3]) 

Z{a) = I dh n{h)[e{a, b) + iri{a, h)]-f{b)z{h) . (15) 



Now Eqs. 



become a system of integral equa- 



tions; still it is simpler than the original equation ([T]) as 
at each value of the continuous parameter a we have only 
three real time-dependent variables. 

Certainly, one can also perform both thermodynamic 
limits simultaneously. Then the ensemble is described by 
the system (jl0lllll3ll5|) . 

Next we study an important case when Eqs. (|10lll|) 
decouple. To this end we represent the fraction in 
Eqs. (|12ll3p as a series 



E (-<e<" 

1=2 



(16) 



where complex constants C° depend only on the distri- 
bution of the constants of motion 



Na 



or CP = 



fc=i 



aa{iP)e''^ dip , (17) 



and we used that = due to Eqs. ((6TT4)) . Obviously, 
the governing equations simplify, if Cf = for ^ > 2 and 
all a, and, hence, 7 = 1. Then the force Z does not 
depend on the phase variable C and Eq. (fTU]) decouples 
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from Eq. pT|) . It is easy to see from Eqs. ((T7| that , 
which are in fact Fourier coefficients of the distribution of 
the constants of motion ip, vanish in the thermodynamic 
hmit of type (i), if (t('0) = l/27r. However, if the number 
of oscillators in a subpopulation Na is finite, then even 
for a uniform spreading of ipk, the discrete sum in (fT7|) 
yields |Cf | = 1, arg(Cf) = for / = Na,2Na, . . ., and 
we get 



7a 



1- 



1 



1 



-^*e*(C''+'/'?) 



(18) 



Thus, the deviation of from unity is exponentially 
small in the size of the subpopulation and, therefore, can 
be neglected for large Na- This is exactly the case, where 
the complex bunch amplitude p is equal to the mean field 
amplitude r, because in (fT2|) 7 = 1. 

Hence, for the uniform distribution of constants of mo- 
tion t/j, ensemble ([1|) admits a simplified description via 
Eq. (jlOp . supplemented by an equation for Za, cither in a 
discrete or in a continuous (for the thermodynamic limit 
of type (ii)) form: 



•"Zb , 



(19) 



Z(a)= dh n{b)e{a,h)-'°''^''-^^ z{b) . (20) 



A relation between the distribution of the original phases 

and the uniform distribution of constants of motion 
follows from Eq. ([5]): one can see that different distribu- 
tions of the phases <j)k, parameterized by different values 
of p, correspond to the uniformly distributed constants 

As a first application of our framework, we apply 
Eqs. ()10lllll2ll5p to the classical Kuramoto problem 
(cf. Q). We set e(a, 6) = e = const, a = 0, use the 
frequency as the subpopulation index, a = lj^ and per- 
form the thermodynamic limit (ii). As a result, in the 
case when 7 = 1 and the variable C, (as well as the 
constants of motion) does not influence the dynamics, 
we obtain exactly Eqs. (|10I20|) . derived recently by Ott 
and Antonsen (OA) 3] under an assumption of a cer- 
tain parameterization of the phase distribution. Consid- 
ering the Lorentzian distribution of natural frequencies 
n{uj) = [tt^oj'^ + 1)]^^ and using analytic properties of 
z{uj) as a function of complex frequency w, OA have cal- 
culated the integral in Eq. ([20|) by the residue of the pole 
at Lo = i and have obtained Z = ez{i). Substituting Z 
into Eq. (|10p ior uj — i, OA derived a closed equation for 
Z = Z/e, i.e. for the usual Kuramoto mean field of the 
whole population: 



i-l + -)Z 



-\Z\^Z 

2' ' 



(21) 



solved it, and in this way obtained explicitely the evolu- 
tion of the mean field. 

From our derivation of the equations of motion we con- 
clude, that the particular ansatz used in ^ corresponds 



to the case of uniformly distributed constants of motion 
"tpk, what is equivalent to vanishing Fourier coefficients 
C/. Next we discuss, what changes if the distribution of 
constants ipk is not uniform, i.e. C/ ^ 0. Let us treat the 
effect of non-vanishing coefficients C; perturbatively, as- 
suming that in the first approximation the OA ansatz is 
valid. Considering for simplicity the effect of C2 7^ only, 
we obtain a correction to the mean field by substituting 
(HI into mi 



AZ : 



^^ z*(c.)(|z(c.)p-l)e'^CMC2(c.) 

7r(a;2 _^ x) ■ ( ) 



Calculation of this integral by the residue yields 



AZ 



ez 



X^)(|z(^)|^-l)e^^^«C2(^) 



(23) 



From Eq. (jll[) it follows that in the first approximation 
C(^) = Co + Therefore AZ oc e^^*. We conclude 
that the contribution from a nonuniform distribution of 
constants V'fc results in an exponentially decaying correc- 
tion to the mean field. The characteristic time scale of 
this decay is 1/2, to be compared with the characteris- 
tic time scale of the evolution of the mean field, which, 
according to (pij) . is (e/2 — 1)"^. Thus, close to criti- 
cality = 2, the approximation of vanishing constants 
Ci works well after short transients; this is not surpris- 
ing as near a bifurcation point the dynamics is typically 
effectively low-dimensional, dominated by a few normal 
modes. Far from criticality the time scale separation is 
not valid and the dynamics is generally high-dimensional. 

As a second example we extend recent results of 
Abrams et al Q. They studied two coupled subpop- 
ulations of identical oscillators, i.e. model H]) with 
a = 1,2, uji = ^^!2^ Ni — N2 and heterogeneous cou- 
pling Ei^i = £2,1 = 2yLt, £1,2 = £2,1 = 2i^, and aa,b = a, 
where 1/ = 1 — fi. Using the OA ansatz [3[ , Abrams et al 
derived equations for the complex order parameters zi.2 
and analyzed the so-called chimera state, where, e.g., the 
first subpopulation is fully synchronized, pi = 1, whereas 
the other one is only partially synchronized, p2 < 1; 
they have found both static and time-periodic solutions 
for p2. With our approach we describe the system ex- 
actly, by writing six Eqs. ([7][5]) for both subpopulations. 
Since we are interested in the chimera state in the second 
subpopulation, the first, synchronous one, is described 
by its phase $1 only. In this case Zi = /ze'*-*^""-* -i- 



,^^e*(-i>2-K/3-a)^ where 



^(P2,1'2)e''''''^^*^' 



-^ ^-^ (24) 



^ ei* 



fe=l 



in. 



Next, we note that the dynamics depends only on the 
phase difference 5 = <I>i — $2 and, hence, write a closed 



4 




FIG. 1: (Color online) Simulation of ensemble (UJ for A'^i — 
N2 = 64, a = 7r/2 — 0.1, and different distributions of con- 
stants of motion tpf. . Mean fields X2, Y2 are defined by 
Eq. Q. (a): fi — 0.6; uniform distribution of ip^^ results 
in the steady state (black plus), cf. whereas nonuniform 
distributions with g = 0.9 and g = 0.7 yield limit cycle so- 
lutions (bold red and blue solid lines, respectively), (b)-(d): 

(2) 

= 0.65; uniform distribution of ^pf. yields a limit cycle solu- 
tion (b), cf. whereas for q = 0.9 (c) and for q = 0.7 (d) we 
observe a new type of the chimera state with a quasiperiodic 
dynamics. 



system of three equations: 

^ = L:J^{^Acos{p-a) + iycosi5-a)) , (25) 

— = —filsma H Asm(p — ajj 

at 2p2 

+ v{A sin(/3 -a-S)- sm{S - a)) , (26) 



2p2 

— „ ^'^ [fj, A sin( (3 — a) -I- i/siniS — a)] 
at 2p2 



(27) 



Following Abrams et al we take a thermodynamic 
Umit N2 — > 00. Next, we take a uniform distribution 
of the constants V'*'^'' (which enter only via the relations 



((24)) ) - we remind that this choice corresponds to the 
restriction imposed by OA on the phase distribution in 
their ansatz. Then from Eq. (j24p it follows A = p2 = r2, 
(3 — and equations for p2 and 6 decouple from Eq. (|27p . 
The obtained Eqs. (|25l26p constitute exactly the system 
analyzed in [5] . For a nonuniform distribution oitpl- , we 
have to analyze the full three-dimensional system (I25II27|) . 
which certainly can exhibit more complex solutions. 

To verify our theoretical prediction, we have performed 
numerical simulations of the ensemble llj for the same 
parameters, where Abrams et al obtained stationary and 
time-periodic solutions, but for different distributions of 
the constants ipk- Namely, we took tpk, uniformly dis- 
tributed in the range —qn < -0^ < qir, where 9 < 1 is 
a parameter. For g = 1 we have reproduced the results 
of 5], while for q < 1 the dynamics attains an additional 
time dependence and becomes periodic and quasiperi- 
odic, respectively (see Fig.[T]). 

In conclusion, we have performed the exact reduction 
of the dynamics of hierarchically organized populations 
of coupled oscillators. Due to the partial integrability, 
only three dynamical variables remain relevant for each 
subpopulation, all other are constants of motion. This 
reduced description is independent of the subpopulation 
sizes and holds also in the thermodynamic limit. In an 
important particular case, when the distribution of the 
constants of motion is uniform, the governing equations 
further decouple and simplify. We have demonstrated 
that this case corresponds to the recently found particu- 
lar ansatz of Ott and Antonsen The analysis of full 
equations has allowed us to estimate corrections to this 
particular solution due to non-uniformity of constants. 
Application of our framework to the model by Abrams 
et al revealed existence of novel, quasiperiodically 
breathing chimera states. 

In this Letter we restricted our attention to the sim- 
plest setups in order to demonstrate applicability of the 
theory. However, the method can be in a straightforward 
way extended to the cases of nonlinearly coupled popula- 
tions 0, externally forced ensembles, etc. In these cases 
even a chaotic dynamics of the global variables can be 
expected. The main limitation of the theory is that the 
coupling in Eq. ([T]) has a sine form. 
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